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Abstract — Sum of Squares programming has been used exten- 
sively over the past decade for the stability analysis of nonlinear 
systems but several questions remain unanswered. In this paper, 
we show that exponential stability of a polynomial vector field 
on a bounded set implies the existence of a Lyapunov function 
which is a sum-of-squares of polynomials. In particular, the 
main result states that if a system is exponentially stable on 
a bounded nonempty set, then there exists an SOS Lyapunov 
function which is exponentially decreasing on that bounded set. 
The proof is constructive and uses the Picard iteration. A bound 
on the degree of this converse Lyapunov function is also given. 
This result implies that semidefinite programming can be used 
to answer the question of stability of a polynomial vector field 
with a bound on complexity. 

I. INTRODUCTION 

Computational and numerical algorithms are extensively 
used in control theory. A particular example is semidefinite 
programming conditions for addressing linear control prob- 
lems, which are formulated as Linear Matrix Inequalities 
(LMIs). Using such tools, several questions on the analysis and 
synthesis of linear systems can be formulated and addressed 
effectively. In fact, ever since the 1990s [1], LMIs have had 
a significant impact in the control field, to the point that once 
the solution of a control problem has been formulated as the 
solution to an LMI, it is considered solved. 

When it comes to nonlinear and infinite-dimensional sys- 
tems, the equivalent problems can be formulated as polynomial 
non-negativity constraints under a Lyapunov framework, but 
these are not, at first glance, as easy to solve. Polynomial 
non-negativity is in fact NP hard. It is for this reason that 
several researchers have looked at alternative tests for non- 
negativity, that are polynomial-time complex to test, and which 
imply non-negativity. One such relaxation is the existence 
of a sum of squares decomposition: the ability to optimize 
over the set of positive polynomials using the sum-of-squares 
relaxation has undoubtedly opened up new ways for addressing 
nonlinear control problems, in much the same way Linear 
Matrix Inequalities are used to address analysis questions 
for linear finite-dimensional systems. However, there remain 
several open questions about how these methods can be used 
to search for Lyapunov functions for nonlinear systems. For 
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references on early work on optimization of polynomials, 
see [2], [3], and [4]. For more recent work see [5] and [6]. For 
a recent review paper, see [7]. Today, there exist a number of 
software packages for optimization over positive polynomials, 
e.g. SOSTOOLS [8] and GloptiPoly [9]. 

At the same time, there are still a number of unanswered 
questions regarding the use of sum of squares as a relaxation 
to nonnegativity and its use for the analysis of nonlinear 
systems. Unanswered questions include, for example, a series 
of questions on controller synthesis and the role of duality 
to convexify this problem, as well as estimating regions of 
attraction of equilibria. On the computation and optimization 
side, it is unclear whether multi-core computing could be used 
for computation, as well as how to take advantage of sparsity 
in semidefinite programming. 

In this paper, we do not consider the problem of computing 
sum-of-squares Lyapunov functions. Such work can be found 
in, e.g. [4], [10]— [12]. Instead, we concentrate on the properties 
of the converse Lyapunov functions for systems of the form 

*(/)=/(*(*)). 

where / : R" — > W is polynomial. In particular, we address 
the question of whether an exponentially stable nonlinear 
system will have a sum-of-squares Lyapunov function which 
establishes this property. This result adds to our previous 
work [13], where we were able to show that exponential sta- 
bility on a bounded set implies the existence of a exponentially 
decreasing polynomial Lyapunov function on that set. 

Work that is relevant to the one presented here includes 
research on continuity properties, see e.g. [14], [15] and [16] 
and the overview in [17]. Infinitely-differentiable functions 
were explored in the work [18], [19]. Other innovative results 
are found in [20] and [21]. The books [22] and [23] treat 
further converse theorems of Lyapunov. Continuity of Lya- 
punov functions is inherited from continuity of the solution 
map with respect to initial condition. An excellent treatment 
of this problem can be found in the text of Arnol'd [24]. 

Unlike the work in [13], this paper is closely tied to systems 
theory as opposed to approximation theory. Our method is to 
take a well-known form of converse Lyapunov function based 
on the solution map and use the Picard iteration to approximate 
the solution map. The advantage of this approach is that if 
the vector field is polynomial, the Picard iteration will also 
be polynomial. Furthermore, the Picard iteration inductively 
retains almost all the properties of the solution map. The result 
is a new form of iterative converse Lyapunov function, Vjt. This 
function is discussed in Section VI. 

The first practical contribution of this paper is to give a 
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bound on the number of decision variables involved in the 
question of exponential stability of polynomial vector fields 
on bounded sets. This is because SOS functions of bounded 
degree can be parameterized by the set of positive matrices 
of fixed size. Furthermore, we note that the question of 
existence of a Lyapunov function with negative derivative is 
convex. Therefore, if the question of polynomial positivity on 
a bounded set is decidable, we can conclude that the problem 
of exponential stability of polynomial vector fields on that 
set is decidable. The further complexity benefit of using SOS 
Lyapunov functions is discussed in Section VIII. 

The main result of the paper is stated and proven in 
Section VI. Preceding the main result is a series of lemmas 
that are used in the proof of the main theorem. In Subsection V 
we show that the Picard iteration is contractive on a certain 
metric space; and in Subsection V-A we propose a new way 
of extending the Picard iteration. In Section V-B we show that 
the Picard iteration approximately retains the differentiability 
properties of the solution map, before we prove the main 
result. The implications of the main result are then explored 
in Section VIII and Section VII. A detailed example is given 
in Section IX. The paper is concluded in Section X. 

II. Main Result 

Before we begin the technical part of the paper, we give a 
simplified version of the main result. 

Theorem 1: Suppose that / is polynomial of degree q and 
that solutions of x = f(x) satisfy 

ll*( f )ll</q*(o)|| e -*' 

for some A > 0, K > 1 and for any x(0) 6 M, where M is 
a bounded nonempty region of radius r. Then there exist 
a,/3,y>0 and a sum-of-squares polynomial V(x) such that 
for any x € M, 

a|W| 2 <V(*)</3||*|| 2 

w(*)7M<-7lWI 2 . 

Further, the degree of V will be less than 
2q( Nk ~ l ) r where k(L,X,K) is any integer such that 
c(Jt) —if Jo (e TL + K(TL) k yK 2 (TL) k < K, and 

\og2K 2 (TL) k 1 



c(ky 

c(k) 2 < 



21 



-K- 



, -(l+c(k))(K + c(k))<-, 

:{\-{2K 2 )-J) 



KL\og2K 2 

and N(L,X,K) is any integer such that NT > '"l^ an d T < 
jj- for some T and where L is a Lipschitz bound on / on B^i - 

III. Sum-of-Squares 

Sum of squares (SOS) methods have been introduced over 
the past decade to allow for the algorithmic solution of 
problems that frequently arise in systems and control the- 
ory, many of which can be formulated as polynomial non- 
negativity constraints that are however difficult to solve. In 
these methods, non-negativity is relaxed to the existence of a 



SOS decomposition, which can be tested using Semidefinite 
programming. 

Consider, for example, the problem of ensuring that a 
polynomial p(x) <G R[x] satisfies p(x) > 0. This problem arises 
naturally when trying to construct Lyapunov functions for 
the stability analysis of dynamical systems, which is the 
topic of this paper. Since ensuring non-negativity is hard [25] 
many researchers have investigated alternative ways to do this. 
In [26], the existence of a Sum of Squares decomposition was 
used for that purpose, which involves the presentation of other 
polynomials pt(x) such that 



p( x ) = L^w 2 



(i) 



Algorithms for ensuring this have appeared in the 1990's [27] 
but it was not until the turn of the century that this was recog- 
nized as being solvable using Semidefinite Programming [28]. 
In particular, (1) can be shown equivalent to the existence of 
a Q h and a vector of monomials Z{x) of degree less than 
or equal half the degree of p{x), such that 

p(x)=Z(x) T QZ(x) 

In the above representation, the matrix Q is not unique, in fact 
it can be represented as 



(2) 



where Q, satisfy Z(x) J Q(L(x) = 0. The search for A; such that 
Q in (2) is such that Q >: is a Linear Matrix Inequality, which 
can be solved using Semidefinite Programming. Moreover, 
if p(x) has unknown coefficients that enter affinely in the 
representation (1), Semidefinite Programming can be used to 
find values for them so that the resulting polynomial is SOS. 

This latter observation can allow us to search for polynomi- 
als that satisfy SOS conditions: the most important example is 
in the construction of Lyapunov functions, which is the topic of 
this paper. For more details, please see [10], [28]. The question 
that we address in this paper is whether Sum of Squares 
Lyapunov functions always exist for locally exponentially 
stable systems. 

IV. Notation and Background 

The core concept we use in this paper is the Picard iteration. 
We use this to construct an approximation to the solution map 
and then use the approximate solution map to construct the 
Lyapunov function. Construction of the Lyapunov function 
will be discussed in more depth later on. 

Denote the Euclidean ball centered at of radius r by B, . 
Consider an ordinary differential equation of the form 



x(t)=f(x(t)), x(0)=x 0> /(0)=0, 



(3) 



where x £ R" and / satisfies appropriate smoothness properties 
for local existence and uniqueness of solutions. The solution 
map is a function <j> which satisfies 

d 



dt 



and 



0(0,*) 
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A. Lyapunov Stability 

The use of Lyapunov functions to prove stability of ordi- 
nary differential equations is well-established. The following 
theorem illustrates the use of Lyapunov functions. 

Definition 2: We say that the system defined by the equa- 
tions in (3) are exponentially stable on X if there exist y, K > 
such that for any io el, 

IWOIl<iq*>lk-* 

for all t > 0. 

Theorem 3 (Lyapunov)-' Suppose there exist constants 
a,j3,y> and a continuously differentiable function V such 
that the following conditions are satisfied for all x S U C M". 

a\\4 2 <V(x) <P\\x\\ 2 

w(i) r /W < -r\\4 2 

Then we have exponential stability of System (3) on 
{x : {y:V(y)<V(x)}cU}. 

B. Fixed-Point Theorems 

Definition 4: Let X be a metric space. A mapping F : X — > 
X is contractive with coefficient d S [0, 1) if 



\Fx — Fy\\ < d \\x —y\ 



ijel. 



The following is a Fixed-Point Theorem. 

Theorem 5 (Contraction Mapping Principle [29]): Let X 
be a complete metric space and let F : X — >• X be a contraction 
with coefficient d. Then there exists a unique y EX such that 



Fy = y. 



Furthermore, for any xq EX, 



F"x -y 



<d k \\x -y\\. 



To apply these results to the existence of the solution map, 
we use the Picard iteration. 

V. Picard Iteration 

We begin by reviewing the Picard iteration. This is the basic 
mathematical tool we will use to define our approximation to 
the solution map and can be found in many texts, e.g. [30]. 

Definition 6: For given T and r, define the complete metric 
space 

X T , r '■= {q{t) : sup t6 [ 7 -] \\q(t) \\ < r, q is continuous. } (4) 
with norm 

k\\x= su p ll?(0ll- 

te[o,T) 

For a fixed x E B r and q E Xj, r , the Picard Iteration [31], 
is defined as 



(Pq){t)±x+ ff(q(s))ds. 
Jo 



In this paper, we also define the Picard iteration iteration on 
functions z(t,x) as 



(ft)(jc > 0=jc+ f f(z(x,s))ds. 
Jo 



We begin by showing that for any radius r, there exists a 
T such that the Picard iteration is contractive on Xj.2r for any 
x E B r . 

Lemma 7: Given r > 0, let T < min{^, ^} where / has 
Lipschitz factor L on Bi r and Q = sup v£ g 2j f(x). Then P : 
Xr,2r — > Xr,2r and there exists some (j> E X T ^ r such that for 
t E [0, T] and x E B r , 



dt 



Ht)=f(H*)), 0(o)=* 



and for any z E Xj^r, 



(j)-P k z 



<(TL) k U-z\\ 



Proof: We first show that for x E B r , P : X T 2r 
q^X T , 2r , then sup ;G[o r] \\q(t)\\ < 2r and so 



■X T2r . If 



\\Pq\\ x = sup 

te[0,T] 



/(?(*)) 



ds 



<\\4+ I \\.f{q{s))\\ds 
Jo 

<r+ [ sup ||/(y)||rfj 

J0 y£B 2r 

<r+TQ< 2r 

Thus we conclude that Pq E Xj,2r- Furthermore, for qi,q2 G 

\\Pq\ - PqiWx = su p / (/M*))-/fc!(>)))<fa 

te[o,T] J o 

< [ sup \\f(q l ( S ))-f(q 2 (s))\\ds 
J o te[0,T] 

<TL sup \\q 1 (s)-q 2 (s)\\=TL\\qi-q2\\ x 

te[0,T] 

Therefore, by the contraction mapping theorem, the Picard 
iteration converges on [0,T] with convergence rate (TL) k . ■ 



A. Picard Extension Convergence Lemma 

In this section we propose an extension to the Picard iter- 
ation approximation. We divide the interval into subintervals 
on which the Picard iteration is guaranteed to converge. On 
each interval, we apply the Picard iteration using the final 
value of the solution estimate from the previous interval as the 
initial condition, x. For a polynomial vector field, the result is 
a piece-wise polynomial approximation which is guaranteed 
to converge on an arbitrary interval - see Figure 1 for an 
illustration. 

Definition 8: Suppose that the solution map exists on t E 
[0,°°) and ||0(f,Jc)|| < ^IWI f° r an Y x € B r . Suppose that / 
has Lipschitz factor L on B\Kr and is bounded on B^Kr with 
bound Q. Given T < min{2|i, £}, let z = and define 

G k Q (t,x) := (P k z)(t,x) 
and for i > 0, define the functions G, recursively as 

G k +l (t,x):=(P k z)(t,G k (T,x)). 
The G k are Picard iterations P k z(t,x) defined on each 
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Fig. 1, The Solution map <j> and the functions G* for k = 1,2,3,4,5 and 
/= 1,2,3 for the system x(t) = —x(t) 3 . The interval of convergence of the 



sub-interval where we substitute the initial condition x i-> 
G\_i(t,x). Define the concatenation of these G k as 

G k (t,x):=G k (t-iT,x) V te[iT,iT + T] and i = 1, • • • ,00. 

If / is polynomial, then the G k are polynomial for any 
i,k and G k is continuously differentiable in x for any k. The 
following lemma provides several properties for the functions 
G k . 

Lemma 9: Given 8 > 0, suppose that the solution map 
0(f ,x) exists on t G [0, 8] and on x G B r . Further suppose that 
||0(f,a,x)|| <A"||jc|| for any x G fi,-- Suppose that / is Lipschitz 
on B^Kr with factor L and bounded with bound Q. Choose 
T < min{2|i:, ±} and integer N > 8/T. Then let G k and G k 
be defined as above. 

Define the function 



N-l 

m = £ 

1=0 



e TL + K{TLf) K A {TLf. 



Given any k sufficiently large so that c(k) < K, then for any 

x GB r , 



sup 

iG[0,5] 



G k (s,x)-<j>(s,x) < c(ife) ||jc 



(5) 



Proof: Suppose x G B r . By assumption, the conditions of 
Lemma 7 are satisfied using r 1 = 2Kr. Let z(t,x) = 0. Define 
the convergence rate d = TL < 1 . By Lemma 7, 



sup 

s€[0,T] 



G k Q (s,x)-<j)(s,x) 



(P k z)(s,x)-^>(s,x) 



= sup 

se[0,T. 

<d k sup \\(j)(s,x)\\ <Kd k \\x 
se[0,T] 



Thus Equation (5) is satisfied on the interval [0, T], We proceed 
by induction. Define 

Ci{k) = j^{e d +Kd k yK 2 d k . 



and suppose that \\G k — 0|l < c,-_i (k) \\x\\ on interval [iT — 
T.iT]. Then 



sup 

se[iT,iT+T] 



G k (s,x) — <j)(s,x) 



= sup 

se\iT,iT+T] 

= sup 

se\iT,iT+T\ 

< sup 

se[iT,iT+T] 

+ sup 

se[iT,iT+T] 



G k (s-iT,x)-(j>(s,x) 
P k z(s - iT, GU(T,x)) -0(5- iT, (/7») 
J*z(a - /r, G{Li (r,*)) - <l>(s - iT, GU{T,x)) 
0(* - iT, G?_i (7») - 0(5 - iT, (iT,*)) 



We treat these final two terms separately. First note that 



GU{T,x) <U(iT,x) 



0(jT,x)-Gf_ 1 (T,x 



<jq*|| +*_!(*) ||*|| 
<(jr+ Ci _i(fc))||*||. 

Since c;_i(A:) < c(fc) < and x G -S,-, we have ||G*_ 1 (T',x)|| < 
(K + Kci-i(k)) \\x\\< 2Kr. Hence 



sup 

se\iT,iT+T} 



P k z(s - iT,GU{T,x)) -<l>(s- iT,GU{T lX )) 



< sup d k 



se[iT,iT+T] 

<Kd k (K + Ci -i(k))\\x\ 



0{s-iT,GU{T,x)) 



< Kd 



GU(T,x) 



Now,ifxGfi r , ||0(i,Jc)|| <Krand since ||G*L, (T,x) \\ <2Kr 
and / is Lipschitz on BuKr, it is well-known that 



sup 

se[iT,iT+T] 



< sup 

se{iT,iT+T] 



0(5 - iT, G k _ { {T,x)) - 0(5 - iT, (i7») 

GU(T,x) - 0(/7» < /Vi (Jt) ||jc|| 



Ms-iT) 



Combining, we conclude that 

sup G k (s — iT,x) — Q(s,x) 

se[iTJT+T] 

< e TL Ci^(k) \\x\\ +Kd k (K + c i - l (k)) \\x\\ 

= {{e d +Kd k )c,_i (Jt) + K 2 d k ) \\x\\ = e,(jfc) ||jc|| . 

Since c,-(fe) < c(fc), and 5 <NT, by induction, we conclude 
that 

sup G^(s,x) — 0(s,x) < c(&) ||jc|| . 
s€[0,5] 



B. Derivative Inequality Lemma 

In this critical lemma, we show that the Picard iteration 
approximately retains the differentiability properties of the 
solution map. The proof is based on induction, with a key step 
based on an approach in [32] (Proof of Thm 4. 14). This lemma 
is then adapted to the extended Picard iteration introduced in 
the previous section. 

Lemma 10: Suppose that the conditions of Lemma 7 are 
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satisfied. Then for any x £ B r and any k > 0, 



satisfied. Then for any x £ B r 



sup 

'6 [0,71 



^(^z)( f ,x) r /(x)-|(P^)( f ,x) 



<££lMI 



Proof: Begin with the identity for k > 1 

(»,*)=*+ I" f((P k - l z)(s,x))ds 
Jo 

= x + Jj{{P k - l z)(s + t,x))ds. 
Then, by differentiating the right-hand side, we get 

= f((P k - 1 z)(0,x)) 

+ fyfiip'-h^s+t^Y^-'z^s+t^ds 

= f((P k - l zMx))+ fvf((P k - l z)(s,x)) T ^-(P k - l z)(s,x)ds 
Jo as 

= /to + f t Vf((P k - 1 z)(s,x)) T ^-(P k - l z)(s,x)d S , 
Jo OS 



-(P k z)(t,x) =1+ f Vf((P k - l z)(s,x)) T ^-(P k - l z)(s,x)ds. 
: Jo ox 



where 4-f denotes partial differentiation of / with respect to 
its /th variable and 

dx 

Now define for k > 1 , 

y k (t,x) := y x (P k z)(t,x) T f(x) - ^(^z)M- 
For > 2, we have 
y*M := ^- x {P k z)( tl x) T f(x) j t (P k z)(t,x) 



ds 



Vf((P«- l z)(s,x)y -(P k - l z)(s,x)ds 

+ fvf((P k - l z)(s 7 x)) T -^-(P k - 1 z)(s,x)f(x)ds 
Jo ox 

: fvf{{P k -h){s,x)f- 

Jo 

±(P^ z )(s,x)f(x)-^-(P k -\)(s,x) 

OX OS 

= f t Vf((P k - 1 z)(s,x)) T y k - 1 (s,x)ds. 
Jo 

This means that since (P k ~ l z)(t ,x) £B2 r , by induction 

sup||y*(0H <T sup Vf{{P k - l z ){t 7 x)) sup ||y t _iM|| 
[o,r] re[o,r] rG[o.r] 

<TL sup H^.iC^x)!! < (T-Z.)^" 1 ) sup ibiMU 

re[0,r] re[o,r] 

Forfe=l, (Pz)(t,x) =x, soyi(f) =/(x) and sup [0T] ||yi(f)|| < 
L||x||. Thus 

sup || w ( f )||<i-^W. 
fe[o,r] ^ 

■ 

We now adapt this lemma to the extended Picard iteration. 
Lemma 11: Suppose that the conditions of Lemma 9 are 



sup 

re [0,71 



l-G k {t,xff{x)-^-G k {t,x) 



(TL) k 

< [ —L(K + c(k))\\x\\ 



Proof: Recall that 



G k (t,x):=Gi(t-iT,x) V te[iT,iT + T] and i = l,--- ,°o. 

and G k +l (t,x) = P k z{t,G k {T,x)) where z = 0. Then for r e 
[?T,zT + T], 

|-G fc (^^) r /to-^G fc (?^) 

- iT,G k (T,x)) + ^P k (t - iT,G k (T lX )) T f(x) 



< 



(TL) k 



G?(7\x) 



As was shown in the proof of Lemma 9, Gj (T,jc) < (K- 
d{k)) \\x\\. Thus for t £ [;T,iT + r], 

(TL)* 



£.&(t,xff(x) - §-G k (t,x) 



Since the c, are non-decreasing, 



< 



sup 

'£[0,5] 



^G*( f ,*) r /(*)-^G*( f ,*) 



-(* + <*(*)) Ml 



<ffi!(tf + c (*))||,||. 



VI. Main Result - A Converse SOS Lyapunov 
Function 

In this section, we combine the previous results to obtain 
a converse Lyapunov function which is also a sum-of-squares 
polynomial. Specifically, we use a standard form of converse 
Lyapunov function and substitute our extended Picard iteration 
for the solution map. Consider the system 



x(t)=f(x(t)), x(0)=x . 



(6) 



Theorem 12: Suppose that / is polynomial of degree q and 
that system (6) is exponentially stable on M with 

||x(t)||<tf||*(0)||«-^ 

where M is a bounded nonempty region of radius r. Then there 
exist a,P,y>0 and a sum-of-squares polynomial V(x) such 
that for any x £ M, 



«W 2 <vto<^W 2 



(7) 

(8) 



W(x) T f(x)<- 7 \\x\\ 2 . 

Further, the degree of V will be less than 2q( Nk ~ l \ where 
k(L,X,K) is any integer such that c(k) < K, 



c(k) 2 
c{k) z < 



\og2K 2 (TL) 



2X 

2 A 



K- 



, -(l+c(k))(K+c(k))<-, (9) 



KL\og2K 2 



(l-(2K 2 )-J). 



(10) 
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where c(k) is defined as 

c(k)= £ ! (e TL +K(TL) k )'K 2 (TL) k , (11) 

i=0 

and N{L,X,K) is any integer such that iVT > and 7" < 

2^ for some T and where L is a Lipschitz bound on / on B^Kr- 

Proof: Define 8 = an d d = TL. By assumption 

N > y- Next, we note that since stability implies /(0) = 0, / is 
bounded on any B, with bound Q = Lr. Thus for B4&, we have 
the bound g = AKrL. By assumption, r < ^ = ^ = 
Therefore, if is defined as above, the conditions of Lemma 9 
are satisfied. Define G k as in Lemma 9. By Lemma 9, if k 
is defined as above, ||G*(s,;c) — Q(s,x) || < c(k) \\Q(s,x)\\ on 
s G [0, 8} and x G B, . 

We propose the following Lyapunov functions, indexed by 

k. 



V k (x) := / G k (s,x) T G k (s,x)ds 
Jo 



We will show that for any k which satisfies Inequali- 
ties (9), (10) and (11), then if we define V(x) = V k (x), we 
have that V satisfies the Lyapunov Inequalities (7) and (8) and 
has degree less than 2q( Nk ~ l \ The proof is divided into four 
parts: 

Upper and Lower Bounded: To prove that V k is a valid 
Lyapunov function, first consider upper boundedness. If x £ B 
and s £ [0,5]. Then 

2 r i 2 

G k {s 7 x) = <j)(s,x)+ G k (s 7 x) T - <j)(s,x) 

< ||^(j,x)|| 2 + II [c?*(j,x) r - ^(j,*) 

As per Lemma 9, ||G fc (,s,x) — <j>(s,x) \\ < c(k) \\<j>(s,x)\\ < 
Kc(k) ||x||. From stability we have ||0(s,x)|| </T||x||. Hence, 

V k (x)=J^ G k (s,x) 2 ds<8K 2 (l+c(k) 2 ) \\xf. 

Therefore the upper boundedness condition is satisfied for any 
k > with J3 = 8K 2 (l + c(k) 2 ) > 0. 

Next we consider the strict positivity condition. First we 
note 

U(s,x)\\ 2 = G k (s,x)+\<i)(s,x)-G k (s,x) ] " 



< 



G k (s,x) + <j)(s,x)-G k (s,x) 



which implies 



G k (s,x) >||<K,,*)|| 2 - 4>M-g*m 



By Lipschitz continuity of /, ||0(s,jic)|| 2 > e 2Ls \ \x\\ 2 and 
\\G k (s,x)-(j>(s,x)\\ <A"c(Jk)||jc||. Thus 

V k (x) = £\\&(s,x)fds> fl-(l-e- 2LS )-8Kc(k) : 

Therefore for k as defined previously, — e~ 2LS ) — 

8Kc(k) 2 > and so the positivity condition holds for some 
a >0. 



Negativity of the Derivative: Next, we prove the derivative 
condition. Recall 

V k (x):= / G k (s,x) T G k (s,x)ds 
Jo 
rt+5 

G k {s-t 1 x) T G k {s-t,x)ds 

it 

then since VV{x(t)) T f{x(t)) = $jV{x(t)), we have by the 
Leibnitz rule for differentiation of integrals, 



dt 



V k (x(t))= G k (8,x(t)) T G k (8,x(t)) ~ G k (0,x(t)YG k (0 7 x(t)) 



rt+S P) 

- / 2G k (s-t,x(t)) T — G k (s-t,x(t))ds 
Jt 01 

+ p 2G k (s -t,x{t)) T ^G k {s -t,x{t))f{x(t))d S 



G k {8xt)) -won 2 



+ [ 5 2G k (s,x(t))' 1 
Jo 



^G k (s,x(t))mt))-yG k ( S ,x(t)) 



ds 



where recall 4-.f denotes partial differentiation of / with 
respect to its z'th variable. As per Lemma 11, we have 



^-G k (t,x(t)) T m))-^G k (t,x(t)) 



<-(K+c(k))\\x(t)\\ 



and as previously noted \\G k (8,x(t)) \\ 2 < (K 2 e~~ 2X ( s ~''> + 
e(fc) 2 )||x(*)|| 2 Also, ||G*(j,*(f))|| <K(l+c(k))\\x(t)\\. We 
conclude that 

|vA.W0)<(^V 2A6 + cW 2 )|Wf)|| 2 -||x(r)|| 2 



+ 28yK(l+c(k))(K + c(k))\\x(t)\\ 2 
I K 2 e- n5 +c(k) 2 -l+28Ky(l +c(k))(K + c(k))j\\x(t)\\ 2 



Therefore, we have strict negativity of the derivative since 

.d k 
~T 

K\og2K 2 d k 



K 2 e- us +c(k) 2 +28^-(l+c(k))(K + c(k)) 



1 



+ c(k) z + 2- 



2X 



-(l+c(k))(K + c(k)) < 1 



Thus %V k (x(t)) < -Y\\x(t)\\ for some y > 0. 

Sum of Squares: Since / is polynomial and z is trivially 
polynomial, (P k z)(s,x) is a polynomial in x and s. Therefore, 
V k (x) is a polynomial for any k > 0. To show that V k is sum- 
of-squares, we first rewrite the function 



N ,-iT 



T-T 



Vk(x) = £ / G\ (s - iT,x) T G k (s - iT,x 



ds. 



Since G\z is a polynomial in all of its arguments, G k (s — 
iT,x) T G k (s — iT,x) is sum-of-squares. It can therefore be 
represented as Rj(x) T Zi(s) T Zj(s)Rj(x) for some polynomial 
vector Ri and matrix of monomial bases Z,. Then 

N ,iT N 

V k (x) = £>(-*) r /. ZiisfZi^dsRiix) = Y,Ri(x) T MiRi(x) 

;— 1 JiT-T 



7 



Where M; = JL_ T Zi(s) T Zi(s)ds > is a constant matrix. This 
proves that is sum-of-squares since it is a sum of sums-of- 
squares. 

We conclude that V = Vj satisfies the conditions of the 
theorem for any k which satisfies Inequalities (9) and (10). 
Degree Bound: Given a k which satisfies the inequality 
conditions on c{k), we consider the resulting degree of G k , 
and hence, of If / is a polynomial of degree q, and y is 
a polynomial of degree d in x, then Py will be a polynomial 
of degree max{l, dq} in x. Thus since z — 0, the degree of 
P k z will be If N > \, then the degree of G\ will be 

q Nk ~ l . Thus the maximum degree of the Lyapunov function 
is 2q( m ~ l \ 

■ 

In the proof of Theorem 12, the integration interval, 8 
was chosen such that the conditions will always be feasible 
for some k > 0. However, this choice may not be optimal. 
Numerical experimentation has shown us that a better degree 
bound may be obtained by varying this parameter in the proof. 
However, the given value is one which we have found to work 
well in the vast majority of cases. 

We conclude this section by commenting on the form of the 
converse Lyapunov function, 

V k (x):= / G k (s,x) T G k (s,x)ds. 
Jo 

Our Lyapunov function is defined using an approximation of 
the solution map. A dual approach to solution of the Hamilton- 
Jacobi-Bellmand Equation was taken in [33] using occupation 
measures instead of Picard iteration. Indeed, the dual space of 
the Sum of Squares Lyapunov functions can be understood in 
terms of moments of such occupation measures [34]. 

As a final note, the proof of Theorem 12 also holds for time- 
varying systems. Indeed the original proof was for this case. 
However, because Sum-of-Squares is rarely used for time- 
varying systems, the result has been simplified to improve 
clarity of presentation. 

A. Numerical Illustration 

To illustrate the degree bound and hence the complexity of 
analyzing a nonlinear system, we plot the degree bound versus 
the exponential convergence rate of the system. For given 
parameters, this bound is obtained by numerically searching 
for the smallest k which satisfies the conditions of Theorem 12. 
The convergence rate parameter can be viewed as a metric 
for the accuracy of the sum-of-squares approach: suppose 
we have a degree bound as a function of convergence rate, 
d(y). If it is not possible to find a sum-of-squares Lyapunov 
function of degree d(y) proving stability, then we know that 
the convergence rate of the system must be less than y. 

As can be seen, as the convergence rate increases, the 
degree bound decreases super-exponentially, so that at 7= 2.4, 
only a quadratic Lyapunov function is required to prove 
stability. For cases where high accuracy is required, the degree 
bound increases quickly; scaling approximately as ei. To 
reduce the complexity of the problem, in come cases less 
conservative bounds on the degree can be found by considering 




.7 1 2 3 4 5 

Exponential Decay Rate 



Fig. 2. Degree bound vs. Exponential Convergence Rate for K = \.2, r = 
Z,= 1, q = 5. Domains X < .1 and A > .7 are plotted separately for clarity. 

the monomial terms in the vector field. If the complexity is still 
unacceptably high, then one can consider the use of parallel 
computing: unlike single-core processing, parallel computing 
power continues to increase exponentially. For a discussion 
on using parallel computing to solve polynomial optimization 
problems, we refer to [35]. 

VII. Quadratic Lyapunov Functions 

In this section, we briefly explore the implications of our 
result for the existence of quadratic Lyapunov functions prov- 
ing exponential stability of nonlinear systems. Specifically, we 
look at when the theorem predicts the existence of a degree 
bound of 2. We first note that when the vector field is linear, 
then q=l, which implies that 2q( Nk ~ 1 ' 1 = 2 independent of N 
and k. Recall N is the number of Picard iterations, k is the 
number of extensions and q is the degree of the polynomial 
vector field, /. Hence an exponentially stable linear system 
has a quadratic Lyapunov function - which is not surprising. 

Instead we consider the case when q ^ 1. In this case, 
for a quadratic Lyapunov function, we require N = k = I - 
a single Picard iteration and no extensions. By examining 
the proof of Theorem 12, we see that if the conditions of 
the theorem are satisfied with N = k = I then V (x) = x T x is 
a Lyapunov function which establishes exponential stability 
of the system. Since this is perhaps the most commonly 
used form of Lyapunov function, it is worth considering how 
conservative it is when applied to nonlinear systems of the 
form 

*(0 =/(*(*))■ 
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Fig. 3. Required decay rate for a quadratic Lyapunov function vs. Lipschitz 
bound for K = 1.2 



In the following corollary we give sufficient conditions on the 
vector field and decay rate for the Lyapunov function x T x to 
prove exponential stability. 

Corollary 1: Suppose that system (6) is exponentially sta- 
ble with 

IWOH<*|Ko)|| e -*' 

for some X > 0, K > 1 and for any x(Q) S M, where M is a 
bounded nonempty region of radius r. Let L be a Lipschitz 
bound for / on B^Kr- Suppose that there exists some ^ > 
8 > such that 

K 2 e~ 2ls +c{ + 2K8L(l+c 1 )(K + ci) < 1 

and K8L < 1, where c\ = K 2 8L. Let V(x) =x T x. Then for 
any i£M, 

V(x)=Vx T f(x)<-p\\x\\ 2 . 

for some j3 > 0. 

Proof: We reconsider the proof of Theorem 12. This time, 
we set N = k = 1 and T = 8 and determine if there exists 
a 8 = T < 2j which satisfies the upper-boundedness, lower- 
boundedness and derivative conditions. Because V(x) = 8x T x, 
the upper and lower boundedness conditions are immediately 
satisfied. The derivative negativity condition is 

K 2 e- 2ls +c(l) 2 +2K8L(l+c(l))(K + c(l)) < 1 

where c(l) = c\ = K 2 8L. This is satisfied by the statement of 
the theorem. ■ 

Note that neither the size of the region we consider nor 
the degree of the vector field plays any role in determining 
the degree bound. To illustrate the conditions for existence 
of a quadratic Lyapunov function, we plot the required decay 
rate vs. the Lipschitz continuity factor in Figure 3 for K = 
1.2. This plot shows that as the Lipschitz continuity of the 
vector field increases (and the field becomes less smooth), the 
conservatism of using the quadratic Lyapunov function x T x 
increases. 



VIII. Implications for Sum- of- Squares 
Programming 

In this section we consider the implications that the above 
results have on Sum of Squares programming. 

A. Bounding the number of decision variables 

Because the set of continuously differentiable functions is 
an infinite-dimensional vector space, the general problem of 
finding a Lyapunov function is an infinite-dimensional feasi- 
bility problem. However, the set of sum-of-squares Lyapunov 
functions with bounded degree is finite-dimensional. The most 
significant implication of our theorem is a bound on the 
number of variables in the problem of determining stability of 
a nonlinear vector field. The nonlinear stability problem can 
now be expressed as a feasibility problem of the following 
form. 

Theorem 13: For a given X, let 2d be the degree bound as- 
sociated with Theorem 12 and define N = 1 ■ If System (6) 
is exponentially stable on M with decay rate X or greater, the 
following is feasible for some a,j3,y > 0. 

Find: P e S w : 

P>0 

a ||x|| 2 < Z(x) T PZ(x) < /3 \\x\\ 2 for all x e M 
V (Z(x) T PZ(x)) T f(x) < -y\\x\\ 2 for all xEM 

where Z(x) be the vector of monomials in x of degree d or 
less. 

Proof: The proof follows immediately from the fact that 
a polynomial V of degree 2d is SOS if and only if there exists 
a P > such that V(x) = Z{x) T PZ(x). ■ 
Our condition bounds the number of variables in the feasibility 
problem associated with Theorem 13. If M is semialgebraic, 
then the conditions in Theorem 13 can be enforced using sum- 
of-squares and the Positivstellensatz [36]. The complexity of 
solving the optimization problem will depend on the complex- 
ity of the Positivstellensatz test. If positivity on a semialgebraic 
set is decidable, as indicated in [37], this implies the question 
of exponential stability on a bounded set is decidable. 

B. Local Positivity 

Another implication of our result is that it reduces the 
complexity of enforcing the positivity constraint. As discussed 
in Section III, semidefinite programming is used to optimize 
over the cone of sums-of-squares of polynomials. There are 
several different ways the stability conditions can be enforced. 
For example, we have the following theorem. 

Theorem 14: Suppose there exist polynomial V and sum- 
of-squares polynomials 51,527*3 and 54 such that the following 
conditions are satisfied for a,y> 0. 

V{x)-a\\x\\ 2 =s l {x)+g(x)s 2 (x) 
~VV(s) T f(x) - 7 ||x|| 2 = s 3 (x)+g(x)s 4 (x) 

Then we have exponential stability of System (6) on 
{x : {y:V(y)<V(x)}cU}. 
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The complexity of the conditions associated with Theo- 
rem 14 is determined by the four sum-of-squares variables, 
S{. Theorem 14 uses the Positivstellensatz multipliers 52 an d 
S4 to ensure that the Lyapunov function need only be positive 
and decreasing on the region X = {x : g(x) > 0}. However, 
as we now know that the Lyapunov function can be assumed 
SOS, we can eliminate the multiplier reducing complexity 
of the problem. 

Theorem 15: Suppose there exist polynomial V and sum- 
of-squares polynomials si, s% and 53 such that the following 
conditions are satisfied for a,y > 0. 

2 



V(je)-a||je|r 



-V(V(x) + a\\x\\ 2 ) T f(x)-y\\x\\ 



si(x) 
■ S2(x) + g(x)s3(x) 



Then we have exponential stability of System (6) for any x(Q) 
such that {y : V(y) < V(jc(0))} C X where X := {x : g(x) > 0}. 

This simplification reduces the size of the SOS variables by 
25% (from 4 to 3). If the semialgebraic set X is defined using 
several polynomials (e.g. a hypercube), then the reduction in 
the number of variables can approach 50% . SDP solvers are 
typically of complexity 0(n 6 ), where n is the dimension of the 
symmetric matrix variable. In the above example we reduced 
n = AN to « = 3N. Thus this simplification can potentially 
decrease computation by a factor of 82%. 

IX. Numerical Example 

In this section, we use the Van-der-Pol oscillator to illustrate 
how the degree bound influences the accuracy of the stability 
test. The zero equilibrium point of the Van-der-Pol oscillator 
is unstable. In reverse-time, however, this equilibrium is stable 
with a domain of attraction bounded by the well-known 
forward-time limit-cycle. The reverse-time dynamics are as 
follows. 

xi{t) = -x 2 {t) 

x 2 (t) = -H(l -x 1 (f) 2 )x 2 (f)+JCi(f) 

For simplicity, we choose ji = I. On a ball of radius r, the 
Lipschitz constant can be found from L = mp xeBr \\Df{x)\\, 
where ||-|| is the maximum singular value norm. We find a 
Lipschitz constant for the Van-der-Pol oscillator on radius 
r = 1 to be 2.1. Numerical simulations indicate K = 1, as 
illustrated in Figure 4. Given these parameters, the degree 
bound plot is illustrated in Figure 5. Note that the choice 
of K = 1 dramatically improves the degree bound. Numerical 
simulation shows the decay rate to be a relatively constant 
X = .542 throughout the unit ball. This is illustrated in 
Figure 6. This gives us an estimate of the degree bound as 
d = 6. 

To find the converse Lyapunov function associated with this 
degree bound we construct the Picard iteration. 



(Pz)(t,x) =x+ f{0)ds=x. 
Jo 

(P 2 z)(t,x)=x+ f f(Pz(s,x))ds 




Fig. 4. Plot of trajectories of the Van-der-Pol Oscillator. We estimate the 
overshoot parameter as K = 1 



Exponential Decay Rate 



Fig. 5. 
Rate 



Degree Bound for the Van-der-Pol Oscillator as a Function of Decay 




f{x)ds 



-f(x)t 



Fig. 6. A semi-log plot of ||jc|| for three trajectories. We estimate X - 
for the Van-der-Pol oscillator 



.542 
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The converse Lyapunov function is 

V(x) = / (P 2 z(s,x)) T (P 2 z(s,x))ds 



/ (x + f(x)s) T (x + f(x)s)ds 
Jo 

f ' 

Jo 



X 



X 

/to 

T ,5 



SI 





si] 


/to. 


si' 


ds 


X 


s 2 I 


./to. 


8 2 /2I 


X 


8 3 /3I 


m 



ds 



lf8 = T = ± 



|, for the Van-der-Pol Oscillator, we get the 
SOS Lyapunov function. 



192 -V(x) 



X 




'48/ 61' 




X 


./to. 




61 I 




./to. 



A" 
/to 



6.93/ 2.45/ 
2.45/ / 



i 2 


X 




/to. 



6.93x + 2.45/(x) 
2.45x + /(x) 



1 T 



6.93x + 2.45/(x) 
2.45x + /(x) 

= (6.93xi - 2.45x 2 ) 2 + (2.45(x l +x\x 2 ) + 4.48x 2 ) 2 
+ (2.45x-x 2 ) 2 + (xi +x\x 2 + 1.45x 2 ) 2 

As per the previous discussion, we use SOSTOOLS to 
verify that this Lyapunov function proves stability. Note that 
we must show the function is decreasing on the ball of 
radius r = .25, as the Lipschitz bound used in the theorem 
is for the ball of radius B\ r . We are able to verify that the 
Lyapunov function is decreasing on the ball of radius r = .25. 
Some level sets of this Lyapunov function are illustrated in 
Figure 7. Through experimentation, we find that when we 
increase the ball to radius r = 1, the Lyapunov function 
is no longer decreasing. We also found that the quadratic 
Lyapunov function V(x) = x x is not decreasing on the ball of 
radius r = .25. Although we believe that our degree bound is 
somewhat conservative, these results indicate the conservatism 
is not excessive. 

To explore the limits of the SOS approach, for degree bound 
2, 4, 6, 8 and 10, we find the maximum unit ball on which 
we are able to find a sum-of-squares Lyapunov function. We 
then use the largest sublevel set of this Lyapunov function on 
which the trajectories decrease as an estimate for the domain 
of attraction of the system. These level sets are illustrated 
in Figure 8. We see that as the degree bound increases, our 
estimate of the domain of attraction improves. 

X. Conclusion 

In this paper, we have used the Picard iteration to construct 
an approximation to the solution map on arbitrarily long 
intervals. We have used this approximation to prove that 
exponential stability of a polynomial vector field on a bounded 
set implies the existence of a Lyapunov function which is a 
sum-of-squares of polynomials with a bound on the degree. 
This implies that the question of exponential stability on a 



0.8 



0.6 



0.2 - 



-0.2 




-0.6 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 

Fig. 7. Level Sets of the converse Lyapunov function, with Ball of radius 
r=.25 




Fig. 8. Best Invariant Region vs. Degree Bound with Limit Cycle 



bounded set may be decidable. Furthermore, the converse 
Lyapunov function we have used in this paper is relatively easy 
to construct given the vector field and may find applications 
in other areas of control. The main result also holds for time- 
varying systems. 

Recently, there has been interest in using semidefinite 
programming for the analysis on nonlinear systems using sum- 
of-squares. This paper clarifies several questions on the appli- 
cation of this method. We now know that exponential stability 
on a bounded set implies the existence of an SOS Lyapunov 
function and we know how complex this function may be. It 
has been recently shown that globally asymptotically stable 
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vector fields do not always admit sum-of-squares Lyapunov 
functions [38]. Still unresolved is the question of the existence 
of polynomial Lyapunov functions for stability of globally 
exponentially stable vector fields. 
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